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Abstract 

We show that the time dependent single electron, nuclear density matrix of an interacting elec- 
tronic system coupled to nuclear degrees of freedom can be exactly reproduced by that of an 
electronic system with arbitrarily specified electron-electron interactions coupled to the same nu- 
clear degrees of freedom, given the initial density matrix of the interacting system. This formalism 
enables the construction of rigorous time dependent density functional theories to study nonadia- 
batic electronic dynamics. We obtain the Runge-Gross and Van Leeuwen theorems as special cases 
in the adiabatic limit. 
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The correlated dynamics of electrons and nuclei drive fundamental mechanisms of chem- 
ical reaction dynamics. Most of these mechanisms involve nonadiabatic electron nuclear 
dynamics, where the adiabatic or Born-Oppenheimer approximation is no longer valid. 
A central problem in understanding correlated electron-nuclear dynamics is the mutual time 
evolution of the electron and nuclear quantum subsystems. Nonadiabatic energy and infor- 
mation transfer between the two subsystems plays a fundamental role in correlated electron- 
nuclear dynamics. Nonadiabatic processes are controlled by nuclear momentum dependent 
derivative couplings between the electronic and nuclear subsystems. Because of their mo- 
mentum dependence, these couplings also manifest themselves through complex geometric 
phase effects, a well known example of which is the Jahn- Teller effect js, 4|. Nonadiabatic 
couplings between the two subsystems are also generally off diagonal in the nuclear subspace 
and introduce nonlocal, history dependent correlations into the electron-nuclear dynamics. 

Due to the complexity of the time evolution, developing ab initio methods to model cor- 
related electron-nuclear dynamics has proved to be a challenging problem. Time dependent 
density functional theory (TDDFT) is a formally exact, successful method developed to un- 
derstand the dynamicalproperties of interacting electron systems. This theory relies on the 
Runge-Gross theorem [5| and its generalization by Van Leeuwen [6], which relate the exact 
single electron density to the single electron time dependent potential. It is also a compu- 
tationally tractable method, and hence is a promising approach that could be extended to 
model correlated electron-nuclear dynamics. TDDFT has recently been extended to model 
certain open quantum systems j3, 

Direct extensions of TDDFT to model correlated electron-nuclear dynamics involve the 
estimation of only single particle densities and more generally diagonal elements of the sys- 



tem density matrix 
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ll| and are based on extensions of DFT to multicomponent 



systems [1^. Although they guarantee the correct estimation of these quantities, they 
are unable to fully account for the dynamical phase correlations due to the nonadiabatic 
derivative couplings between electronic and nuclear subsystems. These correlations are of 

fundamental importance to nonadiabatic energy transfer mechanisms that occur in many 

I 

chemical reaction processes [ll]. This is because nonadiabatic phase correlations that arise 
during the system's time evolution cannot be fully captured by diagonal density matrix ele- 
ments. It is also difficult to derive essential simplifying semiclassical approximations to the 
multicomponent TDDFT framework that are required for tractable computational studies. 
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We propose a formulation that extends the density functional paradigm to model the time 
evolution of interacting electron nuclear systems starting from an arbitrary initial density 
matrix for the entire system. The theory we present is designed to explicitly and exactly 
capture nonadiabatic effects on correlated electron-nuclear dynamics while simplifying the 
effects of electron-electron interactions. Our theory also retains off diagonal elements of the 
density matrix in the nuclear subspace. This property allows our theory to retain informa- 
tion regarding quantum coherences during the system time evolution. Thus, in contrast to 
the multicomponent DFT approach of Refs.(9)-(11), our theory enables an accurate treat- 
ment of dynamical phenomena induced by nonadiabatic transitions during electron-nuclear 
dynamics. This formalism also allows for the straightforward construction of scmiclassical 
approximations to the nuclear dynamics, and also separably includes the effects of nonadia- 
batic couplings, thus allowing for approximations to model nonadiabatic electronic dynamics 
based on traditional TDDFT. To construct this theory, we map the interacting electronic 
system into a reference system where electron-electron interactions can be specified accord- 
ing to convenience, coupled to nuclear degrees of freedom whose interactions are retained. 
Let i?(q, R) be the combined Hamiltonian of the system, with {R, P} labeling the nuclear 
coordinates and {q} the electronic coordinates: 

i^(q,R) = ^ + i/o(q) + nq;R,t) (1) 

y (q; R, t) is the sum of the time dependent single electron-nuclear potential and the purely 
nuclear- nuclear interaction in the Hamiltonian. The density matrix, p{t), of the system 
evolves according to the Liouville equation: 

ihp{t)=[H,p{t)] (2) 

A key issue is that a theory for nonadiabatic dynamics should be able to evaluate time 
dependent properties of electronic dynamics while simultaneously retaining information re- 
garding dynamics in the nuclear subspace. A natural physical variable that contains this 
information is the single electron, reduced nuclear density matrix which is conjugate to 
the electron-nuclear potential. This matrix is made up of diagonal contributions from the 
electronic subspace, corresponding to the electronic probability density, while the nuclear 
part of the density matrix is fully retained. We derive continuity equations for reduced 
electron-nuclear density and current matrices and relate their time evolution directly to the 
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electron-nuclear potential, y (q; R, t). A suitable and physically transparent method to de- 
rive the continuity equations is obtained from analyzing partial Wigner transforms of the 
density matrix in the nuclear subspace: 

p(R,P;t) = J dzexp[-^](R+||p(t)|R-|) (3) 

The quantum Liouville equation can be written in terms of the partial Wigner transformed 
density matrix: 

zn-p(R,P;t) = {i^(R,P),p(R,P;t)} (4) 

The bracket {..} is the Wigner- Weyl-Moyal bracket . The bracket between two 

arbitrary operators A and B is defined as: 

{A, B} = i(R, P) * B{R, P) - B{R, P) * i(R, P) (5) 

The Moyal product,'*', is defined as a bilinear product of the Wigner transforms of A 



and B: 

-hA 

A(K, P) * B(R, P) = A(R, P) exp 



B(R,P) 



I2i 

A(R,P) = Vp- Vr- Vp- Vr (6) 
The Liouville equation for the partial Wigner transformed density, p = p(R, P; t), is: 

= //o(q),p(t)} + {Mq; R),p} (7) 

Writing 3? = {R, P}, single particle electronic variables are defined as traces over p(JR,; t): 

/(x|3fj,t) = Tr^i {/(x)p(t)} ; / ^ {nj} (8) 



The current and density operators are: 

N 

N 



i=l 

N 



2m 

i=l 



We also define the quantities: 



7,(x|g?; t) = Tr,i | [j(x), Hoiq)]pwi?fi; t) 

tL[v] (gfj, t) = ^ . Vr + M^- R, t), } (11) 



Prom these definitions, the equation of motion for the single electron density and current 
are obtained: 

[dt + iL[v]{U,t)]n{^\^,t) = -V,-J(x|K,i) (12) 

[dt + iL[v](^,t)]M?li,t) = -7.(x|3?,i) 

-— Vxv(x;R,t) *n(x|3fi,t) (13) 
m 

In the adiabatic limit, Eqs.(12) and (13) reduce to the traditional continuity and current 
equations for a system of interacting electrons in an external time dependent potential. 
However, Eq.(12) and Eq.(13) are more general and include the effect of electron-nuclear 
nonadiabatic couplings. This is given by the action of the quantum Liouville operator, iL[v] 
on the single electron distribution functions, {n,j}. If the classical limit is taken for the 
nuclei alone, then the time evolution of the electron density has two contributions, one cor- 
responding to the usual continuity equation in the electronic subspace, and the second from 
the time evolution of a classical nuclear subspace off which electrons can scatter. A satisfac- 
tory reference system would need to reproduce correctly the electronic density n(x|3?, t). We 
show that the single electron density matrix and its corresponding current density matrix 
can be reproduced by several hamiltonians which differ from the original hamiltonian in the 
strength of their electron-electron correlation, through the construction of an appropriate 
single electron-nuclear interaction for each such hamiltonian. 

Statement— Let H{q; R) be a Hamiltonian with a single particle electron-nuclear coupling 
and a electron-electron interaction potential given by the pair {f (x; R, t), Vee(r)}. Given an 
initial value, {n(x|3?, 0), j(x|3?, 0)}, and a pair of solutions, {n(x|3?, t), j(x|9?, t)} to the cur- 
rent and continuity equations, Eqs.(12) and (13), a second Hamiltonian H' with an arbi- 
trarily specified electron-electron interaction, Wee{r), can be constructed to reproduce the 
solutions, {n(x|3?, t), j(x|3?, i)} by modifying its single electron-nuclear potential, w{r\^;t). 
Proof —We present a heuristic and constructive derivation. Assume that there exists a 
Hamiltonian, H' with a pair of interactions, {^(x; R), M^ee(q)} from which the single par- 
ticle density n(x|K;t), and current, j(x|3?;t) can be derived. We will show that for an 
arbitrary two particle electron-electron interaction, W^e, a single particle electron-nuclear 
coupling w can be constructed that reproduces the given electron density. By assumption, 
Eq.(12) and Eq.(13) are satisfied for both the Hamiltonians, if(q; R) and i?'(q;R). For 
notational convenience, we define C(x|3'^; ^) = 'y(x|K; t) — w(x|3?; t) and — % — Iw Then 

5 



by a process of subtraction, we find that tfie single particle current and density satisfy: 

7^,(x|5J;t) = M ^^(x|SR;t) + 1 {c,J(x|J?;t)} (14) 

{C(x|3fi,i),n(x|3fi,i)} = (15) 

Two operators A, B star commute, i.e the equation. A* B — B * A = {) is satisfied when A is 
of the form, A = g((x, t); or A = F{B, *), where F{B, *) is a function generated from Moyal 
star product polynomials in B, i.e. A — ^kgk{^,t)B * B * B . . .k times. In addition to 
these solutions, other functions could exist whose star product commutator with B is zero. 
Thus, the solutions to Eq.(15) can be written in the form: 

C(x|3fJ; t) = C[x, t; n(x|K; t), t)] (16) 

Here, fi(t) contains the set of functions which star commute with n, and which cannot be 
generated from the single particle density n. In analogy with quantum mechanical language, 
we can say that the set of functions, n, Q form a complete set of commuting variables. To 
understand the solutions to Eq.(14), we first study it in the approximation where the classical 
limit is taken for the nuclear variables. In this limit Eq.(14) becomes: 

-nVxC+{C(x|3?;t),j(x|3?;t)| =7l(x|3?;t) (17) 
ml ) PB 

Here the bracket, {..}pb = Vr-Vp — Vp-VrIs the classical Poisson bracket defined over 

the nuclear subspace. Physically, the first term of Eq.(17) is the electronic force due to the 

change in the single particle potential, while the second term is the nonadiabatic coupling 

contributions that are a result of varying the electron-nuclear coupling. These two terms add 

to balance the force contributions due to purely electron-electron correlations. In the absence 

of the nonadiabatic coupling between electrons and nuclei, the equation reduces to the force 

acting on the electrons due to the difference in single particle potentials, which exactly 

compensates for the difference in electron-electron interactions, as in traditional TDDFT. 

The nonadiabatic coupling adds a new contribution corresponding to exchange of energy 

between electronic and nuclear degrees of freedom. We show below that if appropriate 

spatial boundary conditions are specified, Eq.(14) can be solved in the mixed quantum 

classical limit (Eq.(17)) at a given time, t — Q and these solutions can be propagated to 

succeeding timesteps in a iterative fashion. To construct an explicit solution of Eq.(17), we 

define: 

{va, Ua) = {VpJ, VrJ}; a^x,y, z. (18) 
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Eq.(17) can be rewritten as: 



V,[n(x I 3fi; t)C] + J2 K'^^- " ^"^^'a " '^xn] C = lu,v{^ I t) (19) 

a 

We define the operator A and a normalized potential a as: 

iA(x|3fJ;i) = [V [i;„Vii„- -MaVp„-] - 



a(x|3fJ; t) = n(x|3fj; t)C(x|3?; t) (20) 
We also assume that the ratios V „j /n and V ,n/n vanish as x — > ±(X). Eq.(17) becomes: 

V.a + ika = 7l (21) 
This equation can be solved exphcitly as follows: Define the operator, T such that it satisfies: 

Vxf = ifA(x I 3?; t) (22) 
Then f is a path ordered exponential in the electronic subspace given by: 

f(x|3fj;t)= J dy Yl '^exp A(z|K;t)-dz (23) 



=^exp[..] in Eq.(23) is the path ordering exponential operator. The operator f is a path 
integral over all paths, ^ : y — > x, ending in (x, t) in the electronic subspace. Since this is 
a symmetric operator it has an inverse, F"^. As a consequence of the boundary conditions, 
r becomes unity as x — > ±oo. Eq.(20) is rewritten as: 

r-'V^ta = (24) 

This can be solved to obtain: 

<7 = f-i(x|K;i) j^t{y\U;t)t.{y\U;t)-dy (25) 

Thus, a formal solution for Eq.(17) has been constructed. By requantizing the classical 
nuclear variables using the Weyl correspondence, a solution can be generated for the fully 
quantum dynamics, Eq.(14). The solution also has a transparent physical interpretation. 
The kernel, r~^(x)r(y) contains the "history" of the interactions between the electronic 
and nuclear degrees of freedom. In the Born-Oppenheimer limit, the path ordering becomes 



unneccesary and the kernel becomes independent of path. Thus, this kernel builds into the 
effective potential, electronic scattering and nonadiabatic dynamical phase correlations. 
We now describe the timepropagation of the solutions to Eq.(14) using a method similar 



to that developed in Ref. \&\. We require as an initial condition, that the initial density 
matrix of the system is known. Consequently, the quantity, 7(x|3?; 0) and its time derivatives 
at t = can be evaluated from the initial value of the density matrix pw{^', 0). Furthermore 
we make the assumption that all the relevant potentials and distribution functions are Taylor 
expandable around the time t = 0. The Taylor expansions for such functions /(x|3fj;i) are 
defined by: 

oo ^ 

/(x|5R;t) = 5^/(™)(0)— (26) 
^-^ ml 

m=0 

Using this Taylor expansion, Eq.(14) can be rewritten as a system of linear difference equa- 
tions: 

k 



7tJ(0) = E [-V.C^'HO) * n(^-')(0) + {C^^\0),p-H0)} ] (27) 

1=0 ^ 



This equation can be solved iteratively. For k = and k = 1, the system of equations 
become: 

7i°2(0) = -VC(°)(0) * n(°)(0) + ^ {c(°\/°)} (28) 



(0) = 1 VC^') (0) * n W (0) + - VC(°) (0) * r2(i) (0) + I 
m m n 



(29) 



Eq.(27) has one unknown, and all the other quantities are known, given the assumptions 
made. It has a unique solution, upto a function of only the nuclear coordinates whose 
Wigner-Moyal bracket with the current j^^\0) is zero. Furthermore, the second equation 
also has a solution with similar properties, given C''°nO)5 since it also satisfies an equation 
of the same mathematical form. This argument can be inductively extended to solve for 
the k-th time derivatives'^ (0), for an arbitrary value of k. Thus, this method provides a 
constructive time dependent solution to Eq.(14). 

We now consider an important special case. When the electronic dynamics of the system 
is adiabatic, the timescales on which the electronic dynamics occur are much faster than the 
timescales for nuclear motion. As a result, the time derivative of the single electron-nuclear 
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density matrix is dominated by electronic dynamics, and the derivatives w.r.t nuclear degrees 
of freedom can be neglected. This implies that the continuity and current equations become: 

dtn{x\^; t) ^ -Vx/(x|5R; t) (30) 

<9j(x|3fJ; t) ^ -7^x13?; t) - — n(x|K; tW^v (31) 

m 

Eqs.(30) and (31) are of the same form as the equations obtained for a purely electronic 
time evolution under an external time varying potential, v. The Runge-Gross theorem [sl 
has been generalized to relate the time evolution of a electron density and current j} to 
an external potential for systems that satisfy equations of this form [sl. Thus, in the adia- 
batic limit. Van Leeuwen's generalization of the Runge-Gross theorem is obtained from this 
formalism. In conclusion, we have shown that the Runge-Gross theorem can be generalized 
to describe the correlated nonadiabatic dynamics of an interacting electron-nuclear system. 
This theory is very general and allows maximal information about the quantum mechanical 
nuclear variables to be retained. Numerically tractable DFTs can be constructed from the 
explicit solutions, Eq.(25),by the application of various semiclassical reduction schemes on 
the nuclear dynamics. For example, for weakly nonadiabatic mixed quantum classical evo- 
lution, the nonadiabatic kernels in Eq.(25) can be evaluated using a Monte Carlo approach 
with the nonadiabatic derivative operators in Eq. (22) being treated as inducing occassional 
"nonadiabatic transitions" according to a suitable Monte Carlo criterion. The construction 
of such functionals will be discussed in a future work. This formalism demonstrates that the 
density functionals required to correctly approximate electronic nonadiabatic dynamics are 
nonlocal functionals and depend on derivatives in the nuclear position and momentum vari- 
ables. The nonlocality is reminiscent of the memory dependence found in standard TDDFT 
jl6|. The effective potentials include quantum mechanical coherence and geometric phase 
factors that are a consequence of the nonadiabatic coupling between nuclear and electronic 
degrees of freedom. 
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